API

FFT-based power spectrum estimator

Implementation of power spectrum estimator, following https://arxiv.org/abs/1704.02357. Apart from interface choices, differences w.r.t. original nbodykit’s implementation https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py are:

  • real space positions are taken at mesh nodes, instead of 0.5 cell shift (matters only for ell > 0 in global line-of-sight)

  • normalization is computed with density obtained by paintaing data/randoms to mesh, instead of relying on \(\bar{n}_{i}\) column in the catalogs

  • FKP weights are treated as other weights

class pypower.fft_power.BasePowerSpectrumStatistics(edges, modes, power_nonorm, nmodes, wnorm=1.0, shotnoise_nonorm=0.0, power_zero_nonorm=None, power_direct_nonorm=None, corr_direct_nonorm=None, sep_direct=None, attrs=None, mpicomm=None)

Bases: BaseClass

Base template power statistic class. Specific power statistic should extend this class.

We recommend accessing power spectrum measurements through get_power(), or __call__() (accessed through my_power_statistic_instance()).

Initialize BasePowerSpectrumStatistics.

Parameters:

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Power spectrum normalization.

  • shotnoise_nonorm (float, default=0.) – Shot noise, without normalization.

  • power_zero_nonorm (float, array, default=0.) – Value of power spectrum at \(k = 0\).

  • power_direct_nonorm (array, default=0.) – Value of pair-count-based ‘direct’ power spectrum estimation, (e.g. PIP and angular upweights correction, eq. 26 of https://arxiv.org/abs/1912.08803), to be added to power_nonorm.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save() and save_txt()) statistics.

classmethod average(*others, weights=None)

Average input power spectra.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters:

  • add_direct (bool, default=True) – Add pair-count-based ‘direct’ power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part only.

  • Results

  • -------

  • power (array)

property k

Wavenumbers.

property kedges

Wavenumber edges.

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3))  # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2))  # restrict second axis to (0, 0.2)
statistic.select((0, 0.3, 0.01))  # rebin to match step size of 0.01 and restrict to (0, 0.3)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

property with_mpi

Whether to use MPI.

class pypower.fft_power.CatalogFFTPower(data_positions1, data_positions2=None, randoms_positions1=None, randoms_positions2=None, shifted_positions1=None, shifted_positions2=None, data_weights1=None, data_weights2=None, randoms_weights1=None, randoms_weights2=None, shifted_weights1=None, shifted_weights2=None, D1D2_twopoint_weights=None, D1R2_twopoint_weights=None, R1D2_twopoint_weights=None, R1R2_twopoint_weights=None, edges=None, ells=(0, 2, 4), los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype='f8', resampler='tsc', interlacing=2, position_type='xyz', weight_type='auto', weight_attrs=None, direct_engine='corrfunc', direct_selection_attrs=None, direct_edges=None, direct_attrs=None, wnorm=None, shotnoise=None, shotnoise_nonorm=None, mode_oversampling=0, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: MeshFFTPower

Wrapper on MeshFFTPower to estimate power spectrum directly from positions and weights.

Initialize CatalogFFTPower, i.e. estimate power spectrum.

Note

To compute the cross-power of samples 1 and 2, provide data_positions2 (and optionally randoms_positions2, shifted_positions2 for the selection function / shifted random catalogs of population 2). To compute (with the correct shot noise estimate) the auto-power of sample 1, but with 2 weights, provide data_positions1 (but no data_positions2, nor randoms_positions2 and shifted_positions2), data_weights1 and data_weights2; randoms_weights2 and shited_weights2 default to randoms_weights1 and shited_weights1, resp.

Warning

In case line-of-sight is not local, one can provide \(\mu\)-edges. In this case, integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges. For example, with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\). In all other cases, integration is performed between \(\mu = -1.0\) and \(\mu = 1.0\).

Parameters:

  • data_positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • data_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second data catalog. See data_positions1.

  • randoms_positions1 (list, array, default=None) – Optionally, positions of the random catalog representing the first selection function. If no randoms are provided, selection function will be assumed uniform.

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

  • shifted_positions1 (array, default=None) – Optionally, in case of BAO reconstruction, positions of the first shifted catalog.

  • shifted_positions2 (array, default=None) – Optionally, in case of BAO reconstruction, positions of the second shifted catalog.

  • data_weights1 (array of shape (N,), default=None) – Optionally, weights in the first data catalog.

  • data_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second data catalog.

  • randoms_weights1 (array of shape (N,), default=None) – Optionally, weights in the first randoms catalog.

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • shifted_weights1 (array, default=None) – Optionally, weights of the first shifted catalog. See data_weights1.

  • shifted_weights2 (array, default=None) – Optionally, weights of the second shifted catalog. See shifted_weights1.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • los (string, array, default='firstpoint') – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize[. If False and input positions do not fit in the the box size, raise a ValueError.

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights and the mesh.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • interlacing (bool, int, default=2) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2).

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with D1D2_twopoint_weights, D1R2_twopoint_weights, etc. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). The method used to compute the normalization of PIP weights can be specified with the keyword “normalization”: “counter” to normalize each pair by eq. 19 of arXiv:1912.08803. In this case “nalways” specifies the number of bits systematically set to 1 minus the number of bits systematically set to 0 (defaulting to 0). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0, nalways = 1.

  • D1D2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first and second data catalogs. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

  • D1R2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first data catalog and second randoms (shifted if provided) catalog. See D1D2_twopoint_weights.

  • R1D2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between second data catalog and first randoms (shifted if provided) catalog. See D1D2_twopoint_weights.

  • R1R2_twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles between first and second randoms (shifted if provided) catalogs. See D1D2_twopoint_weights.

  • direct_engine (string, default='corrfunc') – Engine for direct power spectrum computation (if bitwise weights or twopoint_weights, or direct_selection_attrs); one of [“kdtree”, “corrfunc”].

  • direct_selection_attrs (dict, default={'theta': (0., 2 / 60.)}) – To select pairs to be counted in the direct power spectrum computation, provide mapping between the quantity (string) and the interval (tuple of floats), e.g. {'rp': (0., 20.)} to select pairs with transverse separation ‘rp’ between 0 and 20, {‘theta’: (0., 20.)}` to select pairs with separation angle ‘theta’ between 0 and 20 degrees. One can additionally provide e.g. ‘counts’: [‘D1D2’, ‘D1R2’] to specify direct counts for which the selection is to be applied. If bitwise weights or twopoint_weights are provided, then this direct power is added to the FFT-based power spectrum; else it is subtracted (for e.g. the \(r_{p}\)-cut or \(\theta\)-cut).

  • direct_edges (array, dict) – To compute direct power spectrum by taking the Bessel transform of pair counts (in configuration space), provide these separation bin edges. May be a dictionary, with keys ‘min’ (minimum \(s\), defaults to 0), ‘max’ (maximum \(s\), defaults to maximum separation given input positions), and ‘step’ (defaults to 1).

  • direct_attrs (dict, default=None) – Optional arguments for BaseDirectCorrEngine (if direct_edges is provided), or for BaseDirectPowerEngine; e.g. nthreads.

  • wnorm (float, default=None) – Power spectrum normalization, to use instead of internal estimate obtained with normalization().

  • shotnoise (float, default=None) – Power spectrum shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing unnormalized_shotnoise() by power spectrum normalization.

  • mode_oversampling (int, default=0) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scattered across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.MeshFFTBase

Bases: BaseClass

Class to override to implement FFT-based computations.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.MeshFFTPower(mesh1, mesh2=None, edges=None, ells=(0, 2, 4), los=None, boxcenter=None, compensations=None, wnorm=None, shotnoise=None, shotnoise_nonorm=None, mode_oversampling=0)

Bases: MeshFFTBase

Class that computes power spectrum from input mesh(es), using global or local line-of-sight, following https://arxiv.org/abs/1704.02357. In effect, this class merges nbodykit’s implementation of the global line-of-sight (periodic) algorithm of: https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py with the local line-of-sight algorithm of: https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py

poles

Estimated power spectrum multipoles.

Type:

PowerSpectrumMultipoles

wedges

Estimated power spectrum wedges (if relevant).

Type:

PowerSpectrumWedges

Initialize MeshFFTPower, i.e. estimate power spectrum.

Warning

In case line-of-sight is not local, one can provide \(\mu\)-edges. In this case, integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges. For example, with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\). In all other cases, integration is performed between \(\mu = -1.0\) and \(\mu = 1.0\).

Parameters:

  • mesh1 (CatalogMesh, RealField, ComplexField) – First mesh. If RealField, assumed to be \(1 + \delta\) or \(\bar{n} (1 + \delta)\). In case of ComplexField, assumed to be the FFT of \(\delta\) (or \(1 + \delta\)), i.e. unit density.

  • mesh2 (CatalogMesh, RealField, ComplexField, default=None) – In case of cross-correlation, second mesh, with same size and physical extent (boxsize and boxcenter) that mesh1.

  • edges (tuple, array, dict, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • boxcenter (float, array, default=None) – Box center; defaults to 0. Used only if provided mesh1 and mesh2 are not CatalogMesh.

  • compensations (list, tuple, string, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic’ (resp. ‘cic-sn’) for cic assignment scheme, with (resp. without) interlacing. In case mesh2 is not None (cross-correlation), provide a list (or tuple) of two such strings (for mesh1 and mesh2, respectively). Used only if provided mesh1 or mesh2 are not CatalogMesh.

  • wnorm (float, default=None) – Power spectrum normalization, to use instead of internal estimate obtained with normalization().

  • shotnoise (float, default=None) – Power spectrum shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation or both mesh1 and mesh2 are pmesh.pm.RealField, and in case of auto-correlation is obtained by dividing unnormalized_shotnoise() of mesh1 by power spectrum normalization.

  • mode_oversampling (int, default=0) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.MetaPowerSpectrumStatistics(name, bases, class_dict)

Bases: BaseMetaClass

Metaclass to return correct power spectrum statistic.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

class pypower.fft_power.PowerSpectrumMultipoles(edges, modes, power_nonorm, nmodes, ells, **kwargs)

Bases: BasePowerSpectrumStatistics

Power spectrum multipoles binned in \(k\).

Initialize PowerSpectrumMultipoles.

Parameters:

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • ells (tuple, list.) – Multipole orders.

  • kwargs (dict) – Other arguments for BasePowerSpectrumStatistics.

classmethod average(*others, weights=None)

Average input power spectra.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

classmethod concatenate_proj(*others)

Concatenate input power spectra, along poles.

Parameters:

others (list of PowerSpectrumMultipoles) – List of power spectrum multipoles to be concatenated.

Returns:

new

Return type:

PowerSpectrumMultipoles

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters:

  • add_direct (bool, default=True) – Add direct power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part if even multipoles, imaginary part if odd multipole.

  • Results

  • -------

  • power (array)

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber = k.

property kedges

Wavenumber edges.

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ndim

Return binning dimensionality.

plot(ax=None, fn=None, kw_save=None, show=False)

Plot power spectrum.

Parameters:

  • ax (matplotlib.axes.Axes, default=None) – Axes where to plot samples. If None, takes current axes.

  • fn (string, default=None) – If not None, file name where to save figure.

  • kw_save (dict, default=None) – Optional arguments for matplotlib.figure.Figure.savefig().

  • show (bool, default=False) – Whether to show figure.

Returns:

ax

Return type:

matplotlib.axes.Axes

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims, ells=None)

Restrict statistic to provided coordinate limits and multipoles ells in place.

For example:

statistic.select((0, 0.3))  # restrict first axis to (0, 0.3)
statistic.select((0, 0.2), ells=(0, 2))  # restrict first axis to (0, 0.2), and select monopole and quadrupole
set_power_direct(power_direct_nonorm=None, **kwargs)

Set direct power spectrum.

Parameters:

  • power_direct_nonorm (array, default=None) – Direct estimation of power spectrum to be stored in power_direct_nonorm.

  • **kwargs (dict) – corr_direct_nonorm and sep_direct can be provided, in which case power_direct_nonorm is computed as the Bessel transform of corr_direct_nonorm.

property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

to_wedges(muedges, ells=None)

Transform poles to wedges, with input \(\mu\)-edges.

Parameters:

  • muedges (array) – \(\mu\)-edges.

  • ells (tuple, list, default=None) – Multipole orders to use in the Legendre expansion. If None, all poles are used.

Returns:

wedges – Power spectrum wedges.

Return type:

PowerSpectrumWedges

property with_mpi

Whether to use MPI.

class pypower.fft_power.PowerSpectrumStatistics(*args, statistic='wedge', **kwargs)

Bases: BaseClass

Entry point to power spectrum statistics.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_power.PowerSpectrumWedges(edges, modes, power_nonorm, nmodes, wnorm=1.0, shotnoise_nonorm=0.0, power_zero_nonorm=None, power_direct_nonorm=None, corr_direct_nonorm=None, sep_direct=None, attrs=None, mpicomm=None)

Bases: BasePowerSpectrumStatistics

Power spectrum binned in \((k, \mu)\).

Initialize BasePowerSpectrumStatistics.

Parameters:

  • edges (tuple of ndim arrays) – Edges used to bin power spectrum measurement.

  • modes (array) – Mean “wavevector” (e.g. \((k, \mu)\)) in each bin.

  • power_nonorm (array) – Power spectrum in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Power spectrum normalization.

  • shotnoise_nonorm (float, default=0.) – Shot noise, without normalization.

  • power_zero_nonorm (float, array, default=0.) – Value of power spectrum at \(k = 0\).

  • power_direct_nonorm (array, default=0.) – Value of pair-count-based ‘direct’ power spectrum estimation, (e.g. PIP and angular upweights correction, eq. 26 of https://arxiv.org/abs/1912.08803), to be added to power_nonorm.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save() and save_txt()) statistics.

classmethod average(*others, weights=None)

Average input power spectra.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=True, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters:

  • add_direct (bool, default=True) – Add pair-count-based ‘direct’ power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part only.

  • Results

  • -------

  • power (array)

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber.

property kedges

Wavenumber edges.

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property mu

Cosine angle to line-of-sight.

property muavg

Mode-weighted average \(\mu\).

property muedges

\(\mu\)-edges.

property ndim

Return binning dimensionality.

plot(ax=None, fn=None, kw_save=None, show=False)

Plot power spectrum.

Parameters:

  • ax (matplotlib.axes.Axes, default=None) – Axes where to plot samples. If None, takes current axes.

  • fn (string, default=None) – If not None, file name where to save figure.

  • kw_save (dict, default=None) – Optional arguments for matplotlib.figure.Figure.savefig().

  • show (bool, default=False) – Whether to show figure.

Returns:

ax

Return type:

matplotlib.axes.Axes

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims)

Restrict statistic to provided coordinate limits in place.

For example:

statistic.select((0, 0.3))  # restrict first axis to (0, 0.3)
statistic.select(None, (0, 0.2))  # restrict second axis to (0, 0.2)
statistic.select((0, 0.3, 0.01))  # rebin to match step size of 0.01 and restrict to (0, 0.3)
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

property with_mpi

Whether to use MPI.

pypower.fft_power.find_unique_edges(x, x0, xmin=0.0, xmax=inf, mpicomm=mpi4py.MPI.COMM_WORLD)

Construct unique edges for distribution of Cartesian distances corresponding to coordinates x. Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py.

Parameters:

  • x (list of ndim arrays) – List of ndim (broadcastable) coordinate arrays.

  • x0 (float) – Fundamental coordinate separation.

  • xmin (float, default=0.) – Minimum separation.

  • xmax (float, default=np.inf) – Maximum separation.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The current MPI communicator.

Returns:

edges – Edges, starting at 0, such that each bin contains a unique value of Cartesian distances.

Return type:

array

pypower.fft_power.get_power_statistic(statistic='wedge')

Return BasePowerSpectrumStatistics subclass corresponding to statistic (either ‘wedge’ or ‘multipole’).

pypower.fft_power.get_real_Ylm(ell, m, modules=None)

Return a function that computes the real spherical harmonic of order (ell, m). Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/convpower/fkp.py.

Note

Faster evaluation will be achieved if sympy and numexpr are available. Else, fallback to numpy and scipy’s functions.

Parameters:

  • ell (int) – The degree of the harmonic.

  • m (int) – The order of the harmonic; abs(m) <= ell.

  • modules (str, default=None) – If ‘sympy’, use sympy + numexpr to speed up calculation. If ‘scipy’, use scipy. If None, defaults to sympy if installed, else scipy.

Returns:

Ylm – A function that takes 3 arguments: (xhat, yhat, zhat) unit-normalized Cartesian coordinates and returns the specified Ylm.

Return type:

callable

References

https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form

pypower.fft_power.normalization(*meshs, uniform=False, resampler='cic', cellsize=10.0, fields=None)

Return DESI-like normalization, summing over mesh cells:

\[A = 1/dV \frac{\alpha_{2} \sum_{i} n_{d,1}^{i} n_{r,2}^{i} + \alpha_{1} \sum_{i} n_{d,2}^{i} n_{r,1}^{i}}{2}\]

\(n_{d,1}^{i}\) and \(n_{r,1}^{i}\) are the first data and randoms density meshes, as obtained by painting data \(w_{d}\) and random weights \(w_{r}\) on the same mesh (of cell volume \(dV\)), using the cic assignment scheme. The sum then runs over the mesh cells. \(\alpha_{1} = \sum_{i} w_{d,1}^{i} / \sum_{i} w_{r,1}^{i}\) and \(\alpha_{2} = \sum_{i} w_{d,2}^{i} / \sum_{i} w_{r,2}^{i}\) where the sum of weights is performed over the catalogs. If no randoms are provided, density is supposed to be uniform and mesh1 and mesh2 are assumed to occupy the same physical volume.

Parameters:

  • meshs (list of CatalogMesh, RealField, ComplexField) – Meshes. If RealField, density is assumed to be uniform, mesh.csum()/np.prod(mesh.pm.BoxSize). If ComplexField, assumed to be the FFT of \(\delta\) (or \(1 + \delta\)), i.e. unit density. In case of autocorrelation, provide the two same meshes (or the mesh, and None).

  • uniform (bool, default=False) – Whether to assume uniform selection function (only revelant when meshes are CatalogMesh).

  • resampler (string, ResampleWindow, default='cic') – Particle-mesh assignment scheme. Choices are [‘ngp’, ‘cic’, ‘tsc’, ‘pcs’].

  • cellsize (array, float) – Physical size of mesh cells used to paint meshes (if instance of CatalogMesh).

Returns:

norm – Normalization.

Return type:

float

pypower.fft_power.normalization_from_nbar(nbar, weights=None, data_weights=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Return BOSS/eBOSS-like normalization, summing over \(\bar{n}_{i}\) and weight columns, i.e.:

\[\alpha \sum_{i=0}^{N} w_{i} \bar{n}_{i}\]
Parameters:

  • nbar (array of shape (N,)) – \(\bar{n}_{i}\) (comoving density) column.

  • weights (array of shape (N,), default=None) – Weights, if any.

  • data_weights (array of shape (N,), default=None) – Data weights, to normalize randoms weights by alpha = sum(data_weights)/sum(weights).

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

Returns:

norm – Normalization.

Return type:

float

pypower.fft_power.project_to_basis(y3d, edges, los=(0, 0, 1), ells=None, antisymmetric=False, exclude_zero=False, mode_oversampling=0)

Project a 3D statistic on to the specified basis. The basis will be one of:

  • 2D \((x, \mu)\) bins: \(\mu\) is the cosine of the angle to the line-of-sight

  • 2D \((x, \ell)\) bins: \(\ell\) is the multipole number, which specifies the Legendre polynomial when weighting different \(\mu\) bins.

Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/fftpower.py.

Notes

In single precision (float32/complex64) nbodykit’s implementation is fairly imprecise due to incorrect binning of \(x\) and \(\mu\) modes. Here we cast mesh coordinates to the maximum precision of input edges, which makes computation much more accurate in single precision.

Notes

We deliberately set to 0 the number of modes beyond Nyquist, as it is unclear whether to count Nyquist as \(\mu\) or \(-\mu\) (it should probably be half weight for both). Our safe choice ensures consistent results between hermitian compressed and their associated uncompressed fields.

Notes

The 2D \((x, \ell)\) bins will be computed only if ells is specified. See return types for further details. For both \(x\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. mode mode falls in bin i if edges[i] <= mode < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\).

Warning

Integration over Legendre polynomials for multipoles is performed between the first and last \(\mu\)-edges, e.g. with \(\mu\)-edges [0.2, 0.4, 0.8], integration is performed between \(\mu = 0.2\) and \(\mu = 0.8\).

Parameters:

  • y3d (RealField or ComplexField) – The 3D array holding the statistic to be projected to the specified basis.

  • edges (list of arrays, (2,)) – List of arrays specifying the edges of the desired \(x\) bins and \(\mu\) bins; assumed sorted.

  • los (array_like, default=(0, 0, 1)) – The line-of-sight direction to use, which \(\mu\) is defined with respect to.

  • ells (tuple of ints, default=None) – If provided, a list of integers specifying multipole numbers to project the 2D \((x, \mu)\) bins on to.

  • antisymmetric (bool, default=False) – If y3d is compressed, whether y3d is hermitian-antisymmetric (in which case negative modes are added with weight -1).

  • exclude_zero (bool, default=0) – Whether to exclude zero mode from the sum.

  • mode_oversampling (int, default=0) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”.

Returns:

  • result (tuple) –

    The 2D binned results; a tuple of (xmean2d, mumean2d, y2d, n2d, zero2d), where:

    • xmean2darray_like, (nx, nmu)

      The mean \(x\) value in each 2D bin

    • mumean2darray_like, (nx, nmu)

      The mean \(\mu\) value in each 2D bin

    • y2darray_like, (nx, nmu)

      The mean y3d value in each 2D bin

    • n2darray_like, (nx, nmu)

      The number of values averaged in each 2D bin

    • zero2darray_like, (nmu, )

      Value of the power spectrum at k = 0

  • result_poles (tuple or None) – The multipole results; if ells supplied it is a tuple of (xmean1d, poles, n1d, poles_zero), where:

    • xmean1darray_like, (nx,)

      The mean \(x\) value in each 1D multipole bin

    • polesarray_like, (nell, nx)

      The mean multipoles value in each 1D bin

    • n1darray_like, (nx,)

      The number of values averaged in each 1D bin

    • poles_zeroarray_like, (nell,)

      Value of the power spectrum at k = 0

pypower.fft_power.unnormalized_shotnoise(*meshs)

Return unnormalized shotnoise, as:

\[\sum_{i=1}^{N_{g}} w_{i,g,1} w_{i,g,2} + \alpha^{2} \sum_{i=1}^{N_{r}} w_{i,r,1} w_{i,r,2}\]

Where the sum runs over data (and optionally) shifted/randoms weights of input meshes.

Painting catalog on mesh

Implementation of methods to paint a catalog on mesh; workhorse is CatalogMesh.

pypower.mesh.ArrayMesh(array, boxsize, type='real', nmesh=None, mpiroot=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Turn numpy array into pmesh.pm.Field.

Parameters:

  • array (array) – Mesh numpy array gathered on mpiroot.

  • boxsize (array, float, default=None) – Physical size of the box along each axis.

  • type (str, default='real') – Type of field, ‘real’, ‘complex’, ‘untransposedcomplex’.

  • nmesh (array, int, default=None) – In case mpiroot is None or complex field, mesh size, i.e. number of mesh nodes along each axis.

  • mpiroot (int, default=0) – MPI rank where input array is gathered. If input array is scattered accross all ranks in C ordering, pass mpiroot = None and specify nmesh.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

Returns:

mesh

Return type:

pmesh.pm.Field

class pypower.mesh.CatalogMesh(data_positions, data_weights=None, randoms_positions=None, randoms_weights=None, shifted_positions=None, shifted_weights=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype='f8', resampler='tsc', interlacing=2, position_type='xyz', copy=False, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: BaseClass

Class to paint catalog of positions and weights to mesh.

Initialize CatalogMesh.

Parameters:

  • data_positions (list, array) – Positions in the data catalog. Typically of shape (3, N) or (N, 3).

  • data_weights (array of shape (N,), default=None) – Optionally, data weights.

  • randoms_positions (list, array) – Positions in the randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_weights (array of shape (N,), default=None) – Randoms weights.

  • shifted_positions (array, default=None) – Optionally, in case of BAO reconstruction, positions of the shifted catalog.

  • shifted_weights (array, default=None) – Optionally, in case of BAO reconstruction, weigths of the shifted catalog.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis.

  • boxsize (array, float, default=None) – Physical size of the box along each axis, defaults to maximum extent taken by all input positions, times boxpad.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize[? If False and input positions do not fit in the the box size, raise a ValueError.

  • boxpad (float, default=2.) – When boxsize is determined from positions, take boxpad times the smallest box enclosing positions as boxsize.

  • dtype (string, dtype, default='f8') – The data type to use for the mesh. Input positions and weights are cast to the corresponding (real) precision.

  • resampler (string, ResampleWindow, default='tsc') – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’].

  • interlacing (bool, int, default=2) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2).

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

  • copy (bool, default=False) – If False, avoids copy of positions and weights if they are of (real) type dtype, mpiroot is None, and position_type is “pos” (for positions). Setting to True is only useful if one wants to modify positions or weights that have been passed as input while keeping those attached to the current mesh instance the same.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

clone(data_positions=None, data_weights=None, randoms_positions=None, randoms_weights=None, shifted_positions=None, shifted_weights=None, position_type='xyz', **kwargs)

Clone current instance, i.e. copy and set new positions and weights. Arguments ‘boxsize’, ‘nmesh’, ‘boxcenter’, ‘dtype’, ‘resampler’, ‘interlacing’, ‘mpicomm’, if None, are overriden by those of the current instance.

property compensation

Return dictionary specifying compensation scheme for particle-mesh resampling.

save(filename)

Save to filename.

to_mesh(field=None, dtype=None, compensate=False)

Paint positions/weights to mesh.

Parameters:

  • field (string, default=None) –

    Field to paint to mesh, one of:

    • ”data”: data positions and weights

    • ”shifted”: shifted positions and weights (available only if shifted positions are provided)

    • ”randoms”: randoms positions and weights

    • ”data-normalized_shifted”: shifted positions and weights, renormalized (by alpha)

      such that their sum is same as data weights

    • ”data-normalized_randoms”: randoms positions and weights, renormalized (by alpha)

      such that their sum is same as data weights

    • ”fkp”: FKP field, i.e. data - alpha * (shifted if provided else randoms)

    • None: defaults to “data” if no shifted/randoms, else “fkp”

  • dtype (string, dtype, default='f8') – The data type of the mesh when painting, to override current dtype.

  • compensate (bool, default=False) – Wether to apply compensation for particle-mesh assignment scheme.

Returns:

out – Mesh, with values in “weights” units (not normalized as density).

Return type:

RealField

property with_mpi

Whether to use MPI.

property with_randoms

Whether randoms positions have been provided.

property with_shifted

Whether “shifted” positions have been provided (e.g. for reconstruction).

Approximate window

Implementation of (approximate) window function estimation and convolution. Typically, the window function will be estimated through CatalogSmoothWindow, and window function matrices using PowerSpectrumSmoothWindowMatrix, following https://arxiv.org/abs/2106.06324.

class pypower.smooth_window.CatalogSmoothWindow(randoms_positions1=None, randoms_positions2=None, randoms_weights1=None, randoms_weights2=None, edges=None, projs=None, power_ref=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype=None, resampler=None, interlacing=None, position_type='xyz', weight_type='auto', weight_attrs=None, direct_engine='corrfunc', direct_selection_attrs=None, direct_edges=None, direct_attrs=None, wnorm=None, shotnoise=None, shotnoise_nonorm=None, mode_oversampling=None, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: MeshFFTPower

Wrapper on MeshFFTPower to estimate window function from input random positions and weigths.

Initialize CatalogSmoothWindow, i.e. estimate power spectrum window.

Note

To compute the cross-window of samples 1 and 2, provide randoms_positions2. To compute (with the correct shot noise estimate) the auto-window of randoms 1, but with 2 weights, provide randoms_positions1, randoms_weights1 and randoms_weights2.

Parameters:

  • randoms_positions1 (list, array, default=None) – Positions in the first randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

  • randoms_weights1 (array of shape (N,), default=None) – Optionally, weights in the first randoms catalog.

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu-edges\) (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’).

  • projs (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells. Namely, for maximum theory and output multipole \(\ell_{\mathrm{max}}\), window function multipoles will be computed at wide-angle order 0 up to \(2 \ell_{\mathrm{max}}\) (maximum order yielded by the product of any theory and output Legendre polynomial, see e.g. eq. C8 of https://arxiv.org/pdf/2106.06324.pdf). In addition, if chosen line-of-sight is local (either ‘firstpoint’ or ‘endpoint’), odd poles of the window function will be computed at wide-angle order 1, up to \(2 \ell_{\mathrm{max}} + 1\) (maximum non-zero odd pole of wide angle order 1 generated by even poles up to \(\ell_{\mathrm{max}}\) is \(\ell_{\mathrm{max}} + 1\)). Finally, if any of power_ref.ells is odd, all (even and odd) poles will be computed at wide-angle orders 0 and 1, up to \(2 \ell_{\mathrm{max}}\) and \(2 \ell_{\mathrm{max}} + 1\), respectively.

  • power_ref (PowerSpectrumMultipoles, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for projs, los, boxsize, boxcenter, nmesh, interlacing, resampler, wnorm and mode_oversampling if those are None.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector. If None, defaults to line-of-sight used in estimation of power_ref.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxsize (array, float, default=None) – Physical size of the box along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions. If None, defaults to the value used in estimation of power_ref.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize[. If False and input positions do not fit in the the box size, raise a ValueError.

  • dtype (string, dtype, default=None) – The data type to use for input positions and weights and the mesh. If None, defaults to the value used in estimation of power_ref if provided, else ‘f8’.

  • resampler (string, ResampleWindow, default=None) – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’]. If None, defaults to the value used in estimation of power_ref.

  • interlacing (bool, int, default=None) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2). If None, defaults to the value used in estimation of power_ref.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, else “product_individual”.

    If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of weights if the denominator is zero (defaulting to 0).

  • direct_engine (string, default='corrfunc') – Engine for direct window function computation (if bitwise weights or twopoint_weights, or direct_selection_attrs); one of [“kdtree”, “corrfunc”].

  • direct_selection_attrs (dict, default={'theta': (0., 2 / 60.)}) – To select pairs to be counted in the direct window function computation, provide mapping between the quantity (string) and the interval (tuple of floats), e.g. {'rp': (0., 20.)} to select pairs with transverse separation ‘rp’ between 0 and 20, {‘theta’: (0., 20.)}` to select pairs with separation angle ‘theta’ between 0 and 20 degrees. One can additionally provide e.g. ‘counts’: [‘D1D2’, ‘D1R2’] to specify direct counts for which the selection is to be applied. If bitwise weights or twopoint_weights are provided, then this direct power is added to the FFT-based window function; else it is subtracted (for e.g. the \(r_{p}\)-cut or \(\theta\)-cut).

  • direct_edges (array, dict) – To compute direct window function by taking the Bessel transform of pair counts (in configuration space), provide these separation bin edges. May be a dictionary, with keys ‘min’ (minimum \(s\), defaults to 0), ‘max’ (maximum \(s\), defaults to maximum separation given input positions), and ‘step’ (defaults to 1).

  • direct_attrs (dict, default=None) – Optional arguments for BaseDirectCorrEngine (if direct_edges is provided), or for BaseDirectPowerEngine; e.g. nthreads.

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing unnormalized_shotnoise() by window function normalization.

  • mode_oversampling (int, default=None) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”. If None, defaults to the value used in estimation of power_ref.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

property boxsize

Physical box size.

classmethod concatenate_x(*others, **kwargs)

Concatenate poles. Same argument as PowerSpectrumSmoothWindow.concatenate_x().

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.smooth_window.CorrelationFunctionSmoothWindow(sep, corr, projs, wnorm_ref=None)

Bases: BaseClass

Correlation window function multipoles.

Initialize CorrelationFunctionSmoothWindow.

Parameters:

  • modes (array) – Mean separation.

  • corr (array) – Mean correlation.

  • projs (list) – List of Projection instances or (multipole, wide-angle order) tuples.

save(filename)

Save to filename.

select(rp=(0.0, inf))

Restrict window to contribution from transverse separation in input rp range.

property with_mpi

Whether to use MPI.

class pypower.smooth_window.CorrelationFunctionSmoothWindowMatrix(sep, projsin, projsout=None, window=None, sum_wa=True, default_zero=False, attrs=None)

Bases: BaseMatrix

Class computing matrix for window product in configuration space.

projmatrix

Array of shape (len(self.projsout), len(self.projsin), len(self.x)).

Type:

array

Initialize CorrelationFunctionSmoothWindowMatrix.

Parameters:

  • sep (array) – Input (and ouput) separations.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, sum_wa=sum_wa).

  • window (CorrelationFunctionSmoothWindow, PowerSpectrumSmoothWindow) – Window function to convolve power spectrum with. If a PowerSpectrumSmoothWindow instance is provided, it is transformed to configuration space.

  • sum_wa (bool, default=True) – Whether to perform summation over output wide-angle orders. Always set to True except for debugging purposes.

  • default_zero (bool, default=False) – If a given projection is not provided in window function, set to 0. Else an IndexError is raised.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

static propose_out(projsin, sum_wa=True)

Propose output projections given proposed input projections projsin. If sum_wa is True (typically always the case), return projections with Projection.wa_order set to None (all wide-angle orders have been summed).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders. By default, line-of-sight is chosen as that save in attrs (attrs['los_type']). To override, use input kwargs which will be passed to CorrelationFunctionOddWideAngleMatrix.

run()

Set up transform, i.e. compute matrix:

\[W_{\ell,\ell^{\prime}}^{(n,n^{\prime})}(s) = \delta_{n n^{\prime}} \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If sum_wa is True, or output projout.wa_order is None, sum over \(n\) (always the case except for debugging purposes). For example, see q. D5 and D6 of arXiv:1810.05051.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.smooth_window.PowerSpectrumSmoothWindow(edges, modes, power_nonorm, nmodes, projs, wnorm_ref=None, **kwargs)

Bases: BasePowerSpectrumStatistics

Power spectrum window function multipoles.

Initialize PowerSpectrumSmoothWindow.

Parameters:

  • edges (tuple of ndim arrays) – Edges used to bin window function measurement.

  • modes (array) – Mean “wavenumber” (\(k\)) in each bin.

  • power_nonorm (array) – Power spectrum of randoms in each bin, without normalization.

  • nmodes (array) – Number of modes in each bin.

  • projs (list) – List of Projection instances or (multipole, wide-angle order) tuples.

  • wnorm_ref (float, array, default=None) – Normalization of the reference data power spectrum; defaults to wnorm.

  • kwargs (dict) – Other arguments for BasePowerSpectrumStatistics.

classmethod average(*others, weights=None)

Average input window functions.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

classmethod concatenate_proj(*others)

Concatenate input window functions, along projections.

Parameters:

others (list of PowerSpectrumSmoothWindow) – List of window functions to be concatenated.

Returns:

new

Return type:

PowerSpectrumSmoothWindow

classmethod concatenate_x(*others, select='nmodes', frac_nyq=None)

Concatenate input window functions, along k-coordinates. k-edges and k-coordinates are taken from the first provided window; then expanded by those of other provided windows, if they cover a wider range. Eventually, for each bin (low, high), loop through all windows and select that with the largest number of modes in the given bin, if exists (bins are declared equal when exact floating point matching of (low, high)). In case two windows have the same number of modes in the same bin, the first provided one is selected. Therefore, different results may be obtained when changing the order of input windows.

Note

Typically, you will want to input windows with decreasing box size (largest box size first). Direct window function estimation is taken from the first in others.

Parameters:

  • others (list of PowerSpectrumSmoothWindow) – List of window functions to be concatenated.

  • select (string, default='nmodes') – How to select input windows for each k (if several); ‘nmodes’: select window with highest number of modes.

  • frac_nyq (float, tuple, default=None) – Optionally, fraction of Nyquist frequency where to cut input windows (e.g. 0.8). If a float, the same for all input windows; else a tuple or a list of such fraction for each input window.

Returns:

new

Return type:

PowerSpectrumSmoothWindow

classmethod from_power(power, wa_order=0, **kwargs)

Build window function from input PowerSpectrumMultipoless.

Parameters:

  • power (PowerSpectrumMultipoless) – Power spectrum measurement to convert into PowerSpectrumSmoothWindow.

  • wa_order (int, default=0) – Wide-angle order used for input power spectrum measurement.

Returns:

window

Return type:

PowerSpectrumSmoothWindow

get_power(add_direct=True, remove_shotnoise=True, null_zero_mode=False, divide_wnorm=True, complex=True)

Return power spectrum, computed using various options.

Parameters:

  • add_direct (bool, default=True) – Add direct power spectrum measurement.

  • remove_shotnoise (bool, default=True) – Remove estimated shot noise.

  • null_zero_mode (bool, default=True) – Remove power spectrum at \(k = 0\) (if within edges).

  • divide_wnorm (bool, default=True) – Divide by estimated power spectrum normalization.

  • complex (bool, default=True) – Whether (True) to return the complex power spectrum, or (False) return its real part if even multipoles, imaginary part if odd multipole.

  • Results

  • -------

  • power (array)

property k

Wavenumbers.

property kavg

Mode-weighted average wavenumber = k.

property kedges

Wavenumber edges.

modeavg(axis=0, method=None)

Return average of modes for input axis.

Parameters:

  • axis (int, default=0) – Axis.

  • method (str, default=None) – If None, return average separation from modes. If ‘mid’, return bin mid-points.

Returns:

modeavg – 1D array of size shape[axis].

Return type:

array

property ndim

Return binning dimensionality.

property power

Power spectrum, normalized and with shot noise removed.

rebin(factor=1)

Rebin power spectrum estimation in place, by factor(s) factor. Input factors must divide shape.

save(filename)

Save to filename.

save_txt(filename, fmt='%.12e', delimiter=' ', header=None, comments='# ', **kwargs)

Save power spectrum as txt file.

Warning

Attributes are not all saved, hence there is load_txt() method.

Parameters:

  • filename (str) – File name.

  • fmt (str, default='%.12e') – Format for floating types.

  • delimiter (str, default=' ') – String or character separating columns.

  • header (str, list, default=None) – String that will be written at the beginning of the file. If multiple lines, provide a list of one-line strings.

  • comments (str, default=' #') – String that will be prepended to the header string.

  • kwargs (dict) – Arguments for get_power().

select(*xlims, projs=None)

Restrict statistic to provided coordinate limits and projections projs in place.

For example:

statistic.select((0, 0.3))  # restrict first axis to (0, 0.3)
statistic.select((0, 0.2), projs=[(0, 0), (2, 0)])  # restrict first axis to (0, 0.2), and select monopole and quadrupole
property shape

Return shape of binned power spectrum power.

property shotnoise

Normalized shot noise.

slice(*slices)

Slice statistics in place. If slice step is not 1, use rebin(). For example:

statistic.slice(slice(0, 10, 2), slice(0, 6, 3)) # rebin by factor 2 (resp. 3) along axis 0 (resp. 1), up to index 10 (resp. 6)
statistic[:10:2, :6:3] # same as above, but return new instance.
classmethod sum(*others)

Sum input power spectra, weighted by their wnorm.

Warning

Input power spectra have same edges / number of modes for this operation to make sense (no checks performed).

to_real(**kwargs)

Transform power spectrum window function to configuration space.

Parameters:

kwargs (dict) – Arguments for power_to_correlation_window().

Returns:

window

Return type:

CorrelationFunctionSmoothWindow

property with_mpi

Whether to use MPI.

class pypower.smooth_window.PowerSpectrumSmoothWindowMatrix(kout, projsin, projsout=None, weightsout=None, k=None, kin_rebin=1, kin_lim=(0.0001, 1.0), sep=None, window=None, xy=1, q=0, sum_wa=True, default_zero=False, attrs=None)

Bases: BaseMatrix

Class computing matrix for window convolution in Fourier space.

Initialize PowerSpectrumSmoothWindowMatrix.

Parameters:

  • kout (array, PowerSpectrumMultipoles) – Output wavenumbers or power spectrum estimation to take wavenumbers from. Typically, one would use for kout and weightsout the modes PowerSpectrumMultipoles.k and number of modes PowerSpectrumMultipoles.nmodes of the observed power spectrum, such that output k-modes and their weighting of power spectrum and window matrix match.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to the poles of kout if it is a PowerSpectrumMultipoles, else propose_out(projsin, sum_wa=sum_wa).

  • weightsout (array, list, default=None) – Optionally, list of weights to apply when rebinning output “observed” coordinates. Defaults to the number of modes PowerSpectrumMultipoles.nmodes of kout if it is a PowerSpectrumMultipoles, else no weight.

  • k (array, default=None) – Wavenumber for Hankel transforms; must be log-spaced. If None, use sep and xy instead to determine k.

  • kin_rebin (tuple, default=1) – To save some memory, rebin along input k-coordinates by this factor.

  • kin_lim (tuple, default=(1e-4, 1.)) – To save some memory, pre-cut input k-coordinates to these limits.

  • sep (array, default=None) – Separations for Hankel transforms; must be log-spaced. If None, use k and xy instead to determine sep.

  • window (CorrelationFunctionSmoothWindow, PowerSpectrumSmoothWindow) – Window function to convolve power spectrum with. If a PowerSpectrumSmoothWindow instance is provided, it is transformed to configuration space.

  • xy (float, default=1) – If one of k or sep is None, set it following e.g. xy/sep[::-1].

  • q (int, default=0) – Power-law tilt to regularize Hankel transforms.

  • sum_wa (bool, default=True) – Whether to perform summation over output wide-angle orders. Always set to True except for debugging purposes.

  • default_zero (bool, default=False) – If a given projection is not provided in window function, set to 0. Else an IndexError is raised.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

propose_out(sum_wa=True)

Propose output projections given proposed input projections projsin. If sum_wa is True (typically always the case), return projections with Projection.wa_order set to None (all wide-angle orders have been summed).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders in place. By default, line-of-sight is chosen as that save in attrs (attrs['los_type']). To override, use input kwargs which will be passed to PowerSpectrumOddWideAngleMatrix.

run()

Set matrix. Provided arXiv:2106.06324 eq. 2.5:

\[W_{\ell\ell^{\prime}}^{(n)}(k, k^{\prime}) = \frac{2}{\pi} i^{\ell} (-i)^{\ell^{\prime}} \int ds s^{2} j_{\ell}(ks) j_{\ell^{\prime}}(k^{\prime}s) \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

with \(\ell\) corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(k\) to kout and \(k^{\prime}\) to kin. \(n\) is the wide-angle order proj.wa_order. Yet, to avoid bothering with complex values, we only work with the imaginary part of odd power spectra (input and output). In addition, we include the \(dk^{\prime} k^{\prime 2}\) volume element (arXiv:2106.06324 eq. 2.7). Hence we actually implement:

\[W_{\ell,\ell^{\prime}}^{(n)}(k, k^{\prime}) = dk^{\prime} k^{\prime 2} \frac{2}{\pi} (-1)^{\ell/2} (-1)^{\ell^{\prime}/2} \int ds s^{2} j_{\ell}(ks) j_{\ell^{\prime}}(k^{\prime}s) \sum_{L} C_{\ell \ell^{\prime} L} Q_{L}^{(n)}(s)\]

Note that we do not include \(k^{-n}\) as this factor is included in PowerSpectrumOddWideAngleMatrix.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

pypower.smooth_window.power_to_correlation_window(fourier_window, sep=None, k=None, smooth=None)

Compute correlation window function by taking Hankel transforms of input power spectrum window function.

Parameters:

  • fourier_window (PowerSpectrumSmoothWindow) – Power spectrum window function.

  • sep (array, default=None) – Separations \(s\) where to compute Hankel transform; defaults to inverse of fourier_window wavenumbers.

  • k (array, default=None) – Wavenumbers where to interpolate the window function. If provided, \(k\)-space volume element will be computed as \(4 \pi dk k^{2}\). Else, defaults to fourier_window.k and fourier_window.volume.

  • smooth (float, array, default=None) – If not None, if float, radius of Gaussian smoothing. Else, smoothing kernel, should be the same size as used k (see above).

Returns:

window – Correlation window function.

Return type:

CorrelationFunctionSmoothWindow

pypower.smooth_window.weights_trapz(x)

Return weights for trapezoidal integration.

pypower.smooth_window.wigner3j_square(ellout, ellin, prefactor=True)

Return the coefficients corresponding to the product of two Legendre polynomials, corresponding to \(C_{\ell \ell^{\prime} L}\) of e.g. eq. 2.2 of https://arxiv.org/pdf/2106.06324.pdf, with \(\ell\) corresponding to ellout and \(\ell^{\prime}\) to ellin.

Parameters:

  • ellout (int) – Output order.

  • ellin (int) – Input order.

  • prefactor (bool, default=True) – Whether to include prefactor \((2 \ell + 1)/(2 L + 1)\) for window convolution.

Returns:

  • ells (list) – List of mulipole orders \(L\).

  • coeffs (list) – List of corresponding window coefficients.

Accurate window

Implementation of window function estimation, following https://github.com/cosmodesi/GC_derivations, and https://fr.overleaf.com/read/hpgbwqzmtcxn.

class pypower.fft_window.CatalogFFTWindow(randoms_positions1=None, randoms_positions2=None, randoms_weights1=None, randoms_weights2=None, edgesin=None, projsin=None, edges=None, ells=None, power_ref=None, los=None, nmesh=None, boxsize=None, boxcenter=None, cellsize=None, boxpad=2.0, wrap=False, dtype=None, resampler=None, interlacing=None, position_type='xyz', weight_type='auto', weight_attrs=None, wnorm=None, shotnoise=None, shotnoise_nonorm=None, edgesin_type='smooth', mode_oversampling=None, mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD)

Bases: MeshFFTWindow

Wrapper on MeshFFTWindow to estimate window function from input random positions and weigths.

Initialize CatalogFFTWindow, i.e. estimate power spectrum window matrix.

Note

To compute the cross-window of samples 1 and 2, provide randoms_positions2. To compute (with the correct shot noise estimate) the auto-window of randoms 1, but with 2 weights, provide randoms_positions1, randoms_weights1 and randoms_weights2.

Parameters:

  • randoms_positions1 (list, array, default=None) – Positions in the first randoms catalog. Typically of shape (3, N) or (N, 3).

  • randoms_positions2 (list, array, default=None) – Optionally (for cross-correlation), positions in the second randoms catalog. See randoms_positions1.

  • randoms_weights1 (array of shape (N,), default=None) – Optionally, weights in the first randoms catalog.

  • randoms_weights2 (array of shape (N,), default=None) – Optionally (for cross-correlation), weights in the second randoms catalog.

  • edgesin (dict, array, list) – An array of \(k\)-edges which defines the theory \(k\)-binning; corresponding derivatives will be computed using get_correlation_function_tophat_derivative(); or a dictionary of such array for each theory projection. Else a list of derivatives (callable) of theory correlation function w.r.t. each theory basis vector, e.g. each in \(k\)-bin; or a dictionary of such list for each theory projection.

  • projsin (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells.

  • power_ref (PowerSpectrumMultipoles, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for edges, ells, los, boxsize, boxcenter, nmesh, interlacing, resampler, wnorm and mode_oversampling if those are None.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. bins[i] <= x < bins[i+1]. However, last \(\mu\)-bin is inclusive on both ends: bins[-2] <= mu <= bins[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\). If None, defaults to the edges used in estimation of power_ref.

  • ells (list, tuple, default=(0, 2, 4)) – Output multipole orders.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector. If None, defaults to line-of-sight used in estimation of power_ref.

  • nmesh (array, int, default=None) – Mesh size, i.e. number of mesh nodes along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxsize (array, float, default=None) – Physical size of the box along each axis. If None, defaults to the value used in estimation of power_ref.

  • boxcenter (array, float, default=None) – Box center, defaults to center of the Cartesian box enclosing all input positions. If None, defaults to the value used in estimation of power_ref.

  • cellsize (array, float, default=None) – Physical size of mesh cells. If not None, and mesh size nmesh is not None, used to set boxsize as nmesh * cellsize. If nmesh is None, it is set as (the nearest integer(s) to) boxsize / cellsize.

  • boxpad (float, default=2.) – When boxsize is determined from input positions, take boxpad times the smallest box enclosing positions as boxsize.

  • wrap (bool, default=False) – Whether to wrap input positions in [0, boxsize[. If False and input positions do not fit in the the box size, raise a ValueError.

  • dtype (string, dtype, default=None) – The data type to use for input positions and weights and the mesh. If None, defaults to the value used in estimation of power_ref if provided, else ‘f8’.

  • resampler (string, ResampleWindow, default=None) – Resampler used to assign particles to the mesh. Choices are [‘ngp’, ‘cic’, ‘tcs’, ‘pcs’]. If None, defaults to the value used in estimation of power_ref.

  • interlacing (bool, int, default=None) – Whether to use interlacing to reduce aliasing when painting the particles on the mesh. If positive int, the interlacing order (minimum: 2). If None, defaults to the value used in estimation of power_ref.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, else “product_individual”.

    If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of weights if the denominator is zero (defaulting to 0).

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation and in case of auto-correlation is obtained by dividing unnormalized_shotnoise() by window function normalization.

  • edgesin_type (str, default='smooth') – Technique to transpose edgesin to Fourier space. ‘smooth’ uses get_correlation_function_tophat_derivative(); ‘fourier-grid’ paints edgesin on the Fourier mesh, then takes the FFT.

  • mode_oversampling (int, default=None) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”. If None, defaults to the value used in estimation of power_ref.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=mpi.COMM_WORLD) – The MPI communicator.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_window.MeshFFTWindow(mesh1=None, mesh2=None, edgesin=None, projsin=None, power_ref=None, edges=None, ells=None, los=None, periodic=False, boxcenter=None, compensations=None, wnorm=None, shotnoise=None, shotnoise_nonorm=None, edgesin_type='smooth', mode_oversampling=None, **kwargs)

Bases: MeshFFTPower

Class that computes window function from input mesh(es), using global or local line-of-sight, see:

poles

Window matrix.

Type:

PowerSpectrumFFTWindowMatrix

Initialize MeshFFTWindow.

Parameters:

  • mesh1 (CatalogMesh, RealField, default=None) – First mesh.

  • mesh2 (CatalogMesh, RealField, default=None) – In case of cross-correlation, second mesh, with same size and physical extent (boxsize and boxcenter) that mesh1.

  • edgesin (dict, array, list) – An array of \(k\)-edges which defines the theory \(k\)-binning; corresponding derivatives will be computed (see edgesin_type); or a dictionary of such array for each theory projection. Else a list of derivatives (callable) of theory correlation function w.r.t. each theory basis vector, e.g. each in \(k\)-bin; or a dictionary of such list for each theory projection. If periodic is True, this should correspond to the derivatives of theory power spectrum (instead of correlation function) w.r.t. each theory basis vector, e.g. each in \(k\) bin.

  • projsin (list, default=None) – List of Projection instances or (multipole, wide-angle order) tuples. If None, and power_ref is provided, the list of projections is set to be able to compute window convolution of theory power spectrum multipoles of orders power_ref.ells.

  • power_ref (CatalogFFTPower, MeshFFTPower, PowerSpectrumWedges, PowerSpectrumMultipoles, default=None) – “Reference” power spectrum estimation, e.g. of the actual data. It is used to set default values for edges, ells, los, boxcenter, compensations and wnorm if those are None.

  • edges (tuple, array, default=None) – If los is local (None), \(k\)-edges for poles. Else, one can also provide \(\mu\)-edges (hence a tuple (kedges, muedges)) for wedges. If kedges is None, defaults to edges containing unique \(k\) (norm) values, see find_unique_edges(). kedges may be a dictionary, with keys ‘min’ (minimum \(k\), defaults to 0), ‘max’ (maximum \(k\), defaults to np.pi/(boxsize/nmesh)), ‘step’ (if not provided find_unique_edges() is used to find unique \(k\) (norm) values between ‘min’ and ‘max’). For both \(k\) and \(\mu\), binning is inclusive on the low end and exclusive on the high end, i.e. edges[i] <= x < edges[i+1]. However, last \(\mu\)-bin is inclusive on both ends: edges[-2] <= mu <= edges[-1]. Therefore, with e.g. \(\mu\)-edges [0.2, 0.4, 1.0], the last \(\mu\)-bin includes modes at \(\mu = 1.0\). Similarly, with \(\mu\)-edges [0.2, 0.4, 0.8], the last \(\mu\)-bin includes modes at \(\mu = 0.8\). If None, defaults to the edges used in estimation of power_ref.

  • ells (list, tuple, default=(0, 2, 4)) – Output multipole orders. If None, defaults to the multipoles used in estimation of power_ref.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’), use local (varying) first point (resp. end point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector. If None, defaults to the line-of-sight used in estimation of power_ref.

  • periodic (bool, default=False) – If True, selection function is assumed uniform, periodic. In this case, mesh1 may be None; in this case nmesh and boxsize default to that of power_ref, else may be set with kwargs.

  • boxcenter (float, array, default=None) – Box center; defaults to 0. Used only if provided mesh1 and mesh2 are not CatalogMesh. If None, defaults to the value used in estimation of power_ref.

  • compensations (list, tuple, string, default=None) – Compensations to apply to mesh to (optionally) correct for particle-mesh assignment scheme; e.g. ‘cic’ (resp. ‘cic-sn’) for cic assignment scheme, with (resp. without) interlacing. In case mesh2 is not None (cross-correlation), provide a list (or tuple) of two such strings (for mesh1 and mesh2, respectively). Used only if provided mesh1 or mesh2 are not CatalogMesh.

  • wnorm (float, default=None) – Window function normalization. If None, defaults to the value used in estimation of power_ref, rescaled to the input random weights — which yields a correct normalization of the window function for the power spectrum estimation power_ref. If power_ref provided, use internal estimate obtained with normalization() — which is wrong (the normalization poles.wnorm can be reset a posteriori using the above recipe).

  • shotnoise (float, default=None) – Window function shot noise, to use instead of internal estimate, which is 0 in case of cross-correlation or both mesh1 and mesh2 are pmesh.pm.RealField, and in case of auto-correlation is obtained by dividing unnormalized_shotnoise() of mesh1 by window function normalization.

  • edgesin_type (str, default='smooth') – Technique to transpose edgesin to Fourier space, relevant only if periodic is False. ‘smooth’ uses get_correlation_function_tophat_derivative(); ‘fourier-grid’ paints edgesin on the Fourier mesh (akin to the periodic case), then takes the FFT.

  • mode_oversampling (int, default=None) – If > 0, artificially increase the resolution of the input mesh by a factor 2 * mode_oversampling + 1. In practice, shift the coordinates of the coordinates of the input grid by np.arange(-mode_oversampling, mode_oversampling + 1) along each of x, y, z axes. This reduces “discrete grid binning effects”. If None, defaults to the value used in estimation of power_ref.

  • kwargs (dict) – Arguments for ParticleMesh in case mesh1 is not provided (as may be the case if periodic is True), typically boxsize, nmesh, mpicomm.

property boxsize

Physical box size.

property dtype

Mesh dtype.

property nmesh

Mesh size.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.fft_window.PowerSpectrumFFTWindowMatrix(matrix, xin, xout, projsin, projsout, nmodes, wnorm=1.0, wnorm_ref=None, attrs=None, mpicomm=None)

Bases: BaseMatrix

Window matrix, relating “theory” input to “observed” output.

Initialize PowerSpectrumFFTWindowMatrix.

Parameters:

  • matrix (array) – 2D array representing window matrix.

  • xin (array, list) – List of input “theory” coordinates. If single array, assumed to be the same for all input projections projsin.

  • xout (list) – List of output “theory” coordinates. If single array, assumed to be the same for all output projections projsout.

  • projsin (list) – List of input “theory” projections.

  • projsout (list) – List of output “observed” projections.

  • nmodes (array) – Number of modes in each bin.

  • wnorm (float, default=1.) – Window function normalization.

  • attrs (dict, default=None) – Dictionary of other attributes.

  • mpicomm (MPI communicator, default=None) – The MPI communicator, only used when saving (save()) matrix.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

classmethod from_power(power, xin, projin=(0, 0), **kwargs)

Create window function from input PowerSpectrumMultipoles.

Parameters:
Returns:

matrix

Return type:

PowerSpectrumFFTWindowMatrix

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

resum_input_odd_wide_angle(**kwargs)

Resum odd wide-angle orders. Input kwargs will be passed to PowerSpectrumOddWideAngleMatrix.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

pypower.fft_window.get_correlation_function_tophat_derivative(kedges, ell=0, k=None, **kwargs)

Return a list of callable corresponding to the derivative of the correlation function w.r.t. \(k\)-bins.

Parameters:

  • kedges (array) – \(k\)-edges of the \(k\)-bins.

  • ell (int, default=0) – Multipole order.

  • k (array, default=None) – If None, calculation will be analytic, which will work if ell in [0, 2, 4], or sympy package is installed (such analytic integration with sympy may take several seconds). If not None, this is the \(k\) log-spaced array for numerical FFTlog integration.

  • kwargs (dict) – If k is not None, other arguments for fftlog.PowerToCorrelation.

Returns:

toret – List of callables, taking configuration-space separation s as input.

Return type:

list

Direct power spectrum estimator

Implementation of direct estimation of power spectrum multipoles, i.e. summing over particle pairs. This should be mostly used to sum over pairs at small transverse separations, otherwise the calculation will be prohibitive.

class pypower.direct_power.BaseDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), selection_attrs=None, position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: BaseClass

Direct power spectrum measurement, summing over particle pairs.

Initialize BaseDirectPowerEngine.

Parameters:

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-power, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • weight_attrs – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). The method used to compute the normalization of PIP weights can be specified with the keyword “normalization”: “counter” to normalize each pair by eq. 19 of arXiv:1912.08803. In this case “nalways” specifies the number of bits systematically set to 1 minus the number of bits systematically set to 0 (defaulting to 0). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0, nalways = 1.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

  • selection_attrs (dict, default={'theta': (0., 2 / 60.)}) – To select pairs to be counted, provide mapping between the quantity (string) and the interval (tuple of floats), e.g. {'rp': (0., 20.)} to select pairs with transverse separation ‘rp’ between 0 and 20, {‘theta’: (0., 20.)}` to select pairs with separation angle ‘theta’ between 0 and 20 degrees.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’, ‘midpoint’), use local (varying) first-point (resp. end-point, mid-point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.CorrfuncDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), selection_attrs=None, position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: BaseDirectPowerEngine

Direct power spectrum measurement, using Corrfunc.

Initialize BaseDirectPowerEngine.

Parameters:

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-power, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • weight_attrs – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). The method used to compute the normalization of PIP weights can be specified with the keyword “normalization”: “counter” to normalize each pair by eq. 19 of arXiv:1912.08803. In this case “nalways” specifies the number of bits systematically set to 1 minus the number of bits systematically set to 0 (defaulting to 0). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0, nalways = 1.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

  • selection_attrs (dict, default={'theta': (0., 2 / 60.)}) – To select pairs to be counted, provide mapping between the quantity (string) and the interval (tuple of floats), e.g. {'rp': (0., 20.)} to select pairs with transverse separation ‘rp’ between 0 and 20, {‘theta’: (0., 20.)}` to select pairs with separation angle ‘theta’ between 0 and 20 degrees.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’, ‘midpoint’), use local (varying) first-point (resp. end-point, mid-point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.DirectPower(*args, engine='corrfunc', **kwargs)

Bases: object

Entry point to direct power engines.

Parameters:

  • engine (string, default='kdtree') – Name of direct power engine, one of [‘kdtree’, ‘corrfunc’].

  • args (list) – Arguments for direct power engine, see BaseDirectPowerEngine.

  • kwargs (dict) – Arguments for direct power engine, see BaseDirectPowerEngine.

Returns:

engine

Return type:

BaseDirectPowerEngine

class pypower.direct_power.KDTreeDirectPowerEngine(modes, positions1, positions2=None, weights1=None, weights2=None, ells=(0, 2, 4), selection_attrs=None, position_type='xyz', weight_type='auto', weight_attrs=None, twopoint_weights=None, los='firstpoint', dtype='f8', mpiroot=None, mpicomm=mpi4py.MPI.COMM_WORLD, **kwargs)

Bases: BaseDirectPowerEngine

Direct power spectrum measurement, summing over particle pairs, identified with KDTree.

Initialize BaseDirectPowerEngine.

Parameters:

  • modes (array) – Wavenumbers at which to compute power spectrum.

  • positions1 (list, array) – Positions in the first data catalog. Typically of shape (3, N) or (N, 3).

  • positions2 (list, array, default=None) – Optionally, for cross-power, positions in the second catalog. See positions1.

  • weights1 (array, list, default=None) – Weights of the first catalog. Not required if weight_type is either None or “auto”. See weight_type.

  • weights2 (array, list, default=None) – Optionally, for cross-power, weights in the second catalog. See weights1.

  • ells (list, tuple, default=(0, 2, 4)) – Multipole orders.

  • position_type (string, default='xyz') –

    Type of input positions, one of:

    • ”pos”: Cartesian positions of shape (N, 3)

    • ”xyz”: Cartesian positions of shape (3, N)

    • ”rdd”: RA/Dec in degree, distance of shape (3, N)

    If position_type is “pos”, positions are of (real) type dtype, and mpiroot is None, no internal copy of positions will be made, hence saving some memory.

  • weight_type (string, default='auto') –

    The type of weighting to apply to provided weights. One of:

    • None: no weights are applied.

    • ”product_individual”: each pair is weighted by the product of weights \(w_{1} w_{2}\).

    • ”inverse_bitwise”: each pair is weighted by \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\).

      Multiple bitwise weights can be provided as a list. Individual weights can additionally be provided as float arrays. In case of cross-correlations with floating weights, bitwise weights are automatically turned to IIP weights, i.e. \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1}))\).

    • ”auto”: automatically choose weighting based on input weights1 and weights2,

      i.e. None when weights1 and weights2 are None, “inverse_bitwise” if one of input weights is integer, else “product_individual”.

    In addition, angular upweights can be provided with twopoint_weights. If floating weights are of (real) type dtype and mpiroot is None, no internal copy of weights will be made, hence saving some memory.

  • weight_attrs (dict, default=None) – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). Inverse probability weight is then computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2}))\). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0.

  • weight_attrs – Dictionary of weighting scheme attributes. In case weight_type is “inverse_bitwise”, one can provide “nrealizations”, the total number of realizations (including current one; defaulting to the number of bits in input weights plus one); “noffset”, the offset to be added to the bitwise counts in the denominator (defaulting to 1) and “default_value”, the default value of pairwise weights if the denominator is zero (defaulting to 0). The method used to compute the normalization of PIP weights can be specified with the keyword “normalization”: “counter” to normalize each pair by eq. 19 of arXiv:1912.08803. In this case “nalways” specifies the number of bits systematically set to 1 minus the number of bits systematically set to 0 (defaulting to 0). For example, for the “zero-truncated” estimator (arXiv:1912.08803), one would use noffset = 0, nalways = 1.

  • twopoint_weights (WeightTwoPointEstimator, default=None) – Weights to be applied to each pair of particles. A WeightTwoPointEstimator instance (from pycorr) or any object with arrays sep (separations) and weight (weight at given separation) as attributes (i.e. to be accessed through twopoint_weights.sep, twopoint_weights.weight) or as keys (i.e. twopoint_weights['sep'], twopoint_weights['weight']) or as element (i.e. sep, weight = twopoint_weights).

  • selection_attrs (dict, default={'theta': (0., 2 / 60.)}) – To select pairs to be counted, provide mapping between the quantity (string) and the interval (tuple of floats), e.g. {'rp': (0., 20.)} to select pairs with transverse separation ‘rp’ between 0 and 20, {‘theta’: (0., 20.)}` to select pairs with separation angle ‘theta’ between 0 and 20 degrees.

  • los (string, array, default=None) – If los is ‘firstpoint’ (resp. ‘endpoint’, ‘midpoint’), use local (varying) first-point (resp. end-point, mid-point) line-of-sight. Else, may be ‘x’, ‘y’ or ‘z’, for one of the Cartesian axes. Else, a 3-vector.

  • dtype (string, dtype, default='f8') – The data type to use for input positions and weights.

  • mpiroot (int, default=None) – If None, input positions and weights are assumed to be scatted across all ranks. Else the MPI rank where input positions and weights are gathered.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

run()

Method that computes the power spectrum and set power_nonorm, to be implemented in your new engine.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.direct_power.MetaDirectPower(name, bases, class_dict)

Bases: BaseMetaClass

Metaclass to return correct direct power engine.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

class pypower.direct_power.RegisteredDirectPowerEngine(name, bases, class_dict)

Bases: BaseMetaClass

Metaclass registering BaseDirectPowerEngine-derived classes.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

pypower.direct_power.get_default_nrealizations(weights)

Return default number of realizations given input bitwise weights = the number of bits in input weights plus one.

pypower.direct_power.get_direct_power_engine(engine='corrfunc')

Return BaseDirectPowerEngine-subclass corresponding to input engine name.

Parameters:

engine (string, default='kdtree') – Name of direct power engine, one of [‘kdtree’, ‘corrfunc’].

Returns:

engine – Direct power engine class.

Return type:

type

pypower.direct_power.get_inverse_probability_weight(*weights, noffset=1, nrealizations=None, default_value=0.0, correction=None, dtype='f8')

Return inverse probability weight given input bitwise weights. Inverse probability weight is computed as: \(\mathrm{nrealizations}/(\mathrm{noffset} + \mathrm{popcount}(w_{1} \& w_{2} \& ...))\). If denominator is 0, weight is set to default_value.

Parameters:

  • weights (int arrays) – Bitwise weights.

  • noffset (int, default=1) – The offset to be added to the bitwise counts in the denominator (defaults to 1).

  • nrealizations (int, default=None) – Number of realizations (defaults to the number of bits in input weights plus one).

  • default_value (float, default=0.) – Default weight value, if the denominator is zero (defaults to 0).

  • correction (2D array, default=None) – Optionally, divide weight by correction[nbits1, nbits2] with nbits1, nbits2 the number of non-zero bits in weights.

  • dtype (string, np.dtype) – Type for output weight.

Returns:

weight – IIP weight.

Return type:

array

Wide-angle corrections

Implementation of odd wide-angle matrices: - CorrelationFunctionOddWideAngleMatrix for correlation function - PowerSpectrumOddWideAngleMatrix for power spectrum, following https://arxiv.org/abs/2106.06324.

class pypower.wide_angle.BaseMatrix(value, xin, xout, projsin, projsout, weightsin=None, weightsout=None, vectorout=None, weight=1.0, attrs=None)

Bases: BaseClass

Base class to represent a linear transform of the theory model, from input projections projsin to output projections projsout.

value

2D array representing linear transform. First axis is input, second is output.

Type:

array

xin

List of input “theory” coordinates.

Type:

list

xout

List of output “theory” coordinates.

Type:

list

projsin

List of input “theory” projections.

Type:

list

projsout

List of output “observed” projections.

Type:

list

weightsin

List of weights to apply when rebinning input “theory” coordinates.

Type:

list

weightsout

List of weights to apply when rebinning output “observed” coordinates.

Type:

list

Initialize BaseMatrix.

Parameters:

  • value (array) – 2D array representing linear transform.

  • xin (array, list) – List of input “theory” coordinates. If single array, assumed to be the same for all input projections projsin.

  • xout (list) – List of output “theory” coordinates. If single array, assumed to be the same for all output projections projsout.

  • projsin (list) – List of input “theory” projections.

  • projsout (list) – List of output “observed” projections.

  • weightsin (array, list, default=None) – Optionally, list of weights to apply when rebinning input “theory” coordinates.

  • weightsout (array, list, default=None) – Optionally, list of weights to apply when rebinning output “observed” coordinates.

  • weight (float, default=1.) – Optionally, weight to apply when summing two matrices.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.CorrelationFunctionOddWideAngleMatrix(sep, projsin, projsout=None, wa_orders=1, los='firstpoint', attrs=None)

Bases: BaseMatrix

Class computing matrix for odd wide-angle expansion of the correlation function.

Initialize CorrelationFunctionOddWideAngleMatrix.

Parameters:

  • k (array) – Input (and ouput) separations.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, wa_orders=wa_orders). If output projections have Projection.wa_order None, wide-angle orders are summed over.

  • wa_orders (int, list) – Wide-angle expansion orders. So far order 1 only is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

static propose_out(projsin, wa_orders=1)

Propose output projections (i.e. multipoles at wide-angle order > 0) that can be computed given proposed input projections projsin.

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

run()

Set matrix:

\[M_{\ell\ell^{\prime}}^{(n,n^{\prime})}(s) = - \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right)} \delta_{\ell,\ell - 1} \delta_{n^{\prime},0} + \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right)} \delta_{\ell,\ell + 1} \delta_{n^{\prime},0}\]

if \(\ell\) is odd and \(n = 1\), else:

\[M_{\ell\ell^{\prime}}^{(0,n^{\prime})}(s) = \delta_{\ell,ell^{\prime}} \delta_{n^{\prime},0}\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If output projout.wa_order is None, sum over \(n\) (correct only if no window convolution is accounted for).

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.PowerSpectrumOddWideAngleMatrix(k, projsin, projsout=None, d=1.0, wa_orders=1, los='firstpoint', attrs=None)

Bases: BaseMatrix

Class computing matrix for odd wide-angle expansion of the power spectrum. Adapted from https://github.com/fbeutler/pk_tools/blob/master/wide_angle_tools.py

Initialize PowerSpectrumOddWideAngleMatrix.

Parameters:

  • k (array) – Input (and ouput) wavenumbers.

  • projsin (list) – Input projections.

  • projsout (list, default=None) – Output projections. Defaults to propose_out(projsin, wa_orders=wa_orders). If output projections have Projection.wa_order None, wide-angle orders are summed over.

  • d (float, default=1) – Distance at the effective redshift. Use \(1\) if already included in window functions.

  • wa_orders (int, list) – Wide-angle expansion orders. So far order 1 only is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

  • attrs (dict, default=None) – Dictionary of other attributes.

classmethod average(*others, weights=None)

Average input matrices, weighted by weights or their weight.

Warning

Input matrices have same projsin / projsout / xin / xout to make sense (no checks performed).

classmethod concatenate_proj(*others, axis='in')

Concatenate input matrices along projection axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input projections) or ‘out’ (to stack output projections).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

classmethod concatenate_x(*others, axis='in')

Concatenate input matrices along x-axis axis.

Parameters:

  • others (BaseMatrix) – Matrices to concatenate.

  • axis (string, default='in') – Should be either ‘in’ (to stack input x) or ‘out’ (to stack output x).

Returns:

matrix – New matrix, of same type as others[0].

Return type:

BaseMatrix

dot(array, unpack=False)

Apply linear transform to input array. If unpack is True, return “unpacked” array, i.e. a list of arrays corresponding to projsout.

index_x(axis='out', xlim=None, projs=None, concatenate=True)

Return indices for given input x-limits and projs.

Parameters:

  • axis (str, default='out') – Axis, ‘in’ or ‘out’.

  • xlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projs (list, default=None) – List of input projections to return indices for. Defaults to projsin if axis == 'in', projsout if axis == 'out'.

  • concatenate (bool, default=True) – If False, return list of indices, for each input projection.

Returns:

indices

Return type:

array, list

static join(*others)

Join input matrices, i.e. dot them, optionally selecting input and output projections such that they match.

property nprojs

Number of input, output projections.

property nx

Tuple of list of length of input and output coordinates.

pack(matrix)

Set matrix from “unpacked” matrix, i.e. from a list of lists of matrices, where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)]. See unpacked().

prod_proj(array, axes=('in', 0), projs=None)

Multiply current matrix by input array along input axes, projection-wise, i.e. a same operation is applied for all coordinates of a given (input projection, output projection) block.

Parameters:

  • array (1D or 2D array) – Array to multiply matrix with.

  • axes (string, tuple) – Tuple of axes to sum over (axis in current matrix (“in” or “out”)), axis in input array). If array is 1D, one can just provide the axis in current matrix (“in” or “out”).

propose_out(wa_orders=1)

Propose output projections (i.e. multipoles at wide-angle order > 0) that can be computed given proposed input projections projsin.

rebin_x(factorin=1, factorout=1, projsin=None, projsout=None, statistic=None)

Rebin current instance. Internal weights weightsin, weightsout, if not None, are applied.

Parameters:

  • factorin (int, default=1) – Rebin matrix along input coordinates by this factor.

  • factorout (int, default=1) – Rebin matrix along output coordinates by this factor.

  • projsin (list, default=None) – List of input projections to apply rebinning to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply rebinning to. Defaults to projsout.

  • statistic (string, callable, default=None) – Operation to apply when performing rebinning. Defaults to average along input coordinates and sum along output coordinates.

run()

Set matrix:

\[M_{\ell\ell^{\prime}}^{(n,n^{\prime})}(k) = - \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right) d} \delta_{\ell,\ell - 1} \delta_{n^{\prime},0} \left[\frac{\ell - 1}{k} - \partial_{k} \right] - \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right) d} \delta_{\ell,\ell + 1} \delta_{n^{\prime},0} \left[ \frac{\ell + 2}{k} + \partial_{k} \right]\]

if \(\ell\) is odd and \(n = 1\), else:

\[M_{\ell\ell^{\prime}}^{(0,n^{\prime})}(k) = \delta_{\ell,ell^{\prime}} \delta_{n^{\prime},0}\]

with \(\ell\) multipole order corresponding to projout.ell and \(\ell^{\prime}\) to projin.ell, \(n\) wide angle order corresponding to projout.wa_order and \(n^{\prime}\) to projin.wa_order. If output projout.wa_order is None, sum over \(n\) (correct only if no window convolution is accounted for). Derivatives \(\partial_{k}\) are computed with finite differences, see arXiv:2106.06324 eq. 3.3.

save(filename)

Save to filename.

select_proj(projsin=None, projsout=None, **kwargs)

Restrict current instance to provided projections.

Parameters:

  • projsin (list, default=None) – List of input projections to restrict to. Defaults to projsin. If one projection is not in projsin, add a new column to matrix, setting a diagonal matrix where input and output projection match (if the case); see xin.

  • projsout (list, default=None) – List of output projections to restrict to. Defaults to projsout. If one projection is not in projsout, add a new row to matrix, setting a diagonal matrix where input and output projection match (if the case); see xout.

  • kwargs (dict) – In case a new input/output projection must be added, xin/xout for this projection.

select_x(xinlim=None, xoutlim=None, projsin=None, projsout=None)

Restrict current instance to provided coordinate limits in place.

Parameters:

  • xinlim (tuple, default=None) – Restrict input coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • xoutlim (tuple, default=None) – Restrict output coordinates to these (min, max) limits. Defaults to (-np.inf, np.inf).

  • projsin (list, default=None) – List of input projections to apply limits to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply limits to. Defaults to projsout.

slice_x(slicein=None, sliceout=None, projsin=None, projsout=None)

Slice matrix in place. If slice step is not 1, use rebin().

Parameters:

  • slicein (slice, default=None) – Slicing to apply to input coordinates, defaults to slice(None).

  • sliceout (slice, default=None) – Slicing to apply to output coordinates, defaults to slice(None).

  • projsin (list, default=None) – List of input projections to apply slicing to. Defaults to projsin.

  • projsout (list, default=None) – List of output projections to apply slicing to. Defaults to projsout.

unpacked(axis=None)

Return unpacked matrix, a list of lists of matrices where block for output projection projout and input projection projin is obtained through matrix[self.projsout.index(projout)][self.projsin.index(projin)].

property with_mpi

Whether to use MPI.

class pypower.wide_angle.Projection(ell, wa_order='default', default_wa_order=0)

Bases: BaseClass

Class representing a “projection”, i.e. multipole and wide-angle expansion order.

ell

Multipole order.

Type:

int

wa_order

Wide-angle order.

Type:

int, None

Initialize Projection.

Parameters:

  • ell (int) – Multipole order.

  • wa_order (int, None, default='default') – Wide-angle order. If ‘default’, defaults to default_wa_order.

  • default_wa_order (int, default=0) – Default wide-angle order to use if wa_order is ‘default’.

clone(**kwargs)

Clone current projection, optionally updating ell and wa_order (using kwargs).

latex(inline=False)

Return latex string for current projection. If inline is True, add surrounding dollar $ signs.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

pypower.wide_angle.derivative_matrix_nonuniform(x)

Second-order finite-difference matrix for d/dx on an arbitrary 1D grid.

\(f'(x_i) \approx -\frac{h_{i+1}}{h_i (h_i + h_{i+1})} f(x_{i-1}) + \frac{h_{i+1} - h_i}{h_i h_{i+1}} f(x_i) + \frac{h_i}{h_{i+1} (h_i + h_{i+1})} f(x_{i+1})\) with \(h_i = x_i - x_{i-1}\) and \(h_{i+1} = x_{i+1} - x_i\).

pypower.wide_angle.odd_wide_angle_coefficients(ell, wa_order=1, los='firstpoint')

Compute coefficients of odd wide-angle expansion, i.e.:

\[- \frac{\ell \left(\ell - 1\right)}{2 \ell \left(2 \ell - 1\right)}, \frac{\left(\ell + 1\right) \left(\ell + 2\right)}{2 \ell \left(2 \ell + 3\right)}\]

For the first point line-of-sight. See https://fr.overleaf.com/read/hpgbwqzmtcxn. A minus sign is applied on both factors if los is ‘endpoint’.

Parameters:

  • ell (int) – (Odd) multipole order.

  • wa_order (int, default=1) – Wide-angle expansion order. So far only order 1 is supported.

  • los (string) –

    Choice of line-of-sight, either:

    • ’firstpoint’: the separation vector starts at the end of the line-of-sight

    • ’endpoint’: the separation vector ends at the end of the line-of-sight.

Returns:

  • ells (list) – List of multipole orders of correlation function.

  • coeffs (list) – List of coefficients to apply to correlation function multipoles corresponding to output ells.

Utilities

A few utilities.

class pypower.utils.BaseClass

Bases: object

Base class that implements copy(). To be used throughout this package.

save(filename)

Save to filename.

property with_mpi

Whether to use MPI.

class pypower.utils.BaseMetaClass(name, bases, class_dict)

Bases: type

Metaclass to add logging attributes to BaseClass derived classes.

mro()

Return a type’s method resolution order.

set_logger()

Add attributes for logging:

  • logger

  • methods log_debug, log_info, log_warning, log_error, log_critical

pypower.utils.cartesian_to_sky(positions, wrap=True, degree=True, dtype=None)

Transform cartesian coordinates into distance, RA, Dec.

Parameters:

  • positions (array of shape (3, N), list of 3 arrays) – Positions in cartesian coordinates.

  • wrap (bool, default=True) – Whether to wrap RA in \([0, 2 \pi]\).

  • degree (bool, default=True) – Whether RA, Dec are in degrees (True) or radians (False).

Returns:

rdd – Right ascension, declination and distance.

Return type:

list of 3 arrays

pypower.utils.distance(positions)

Return cartesian distance, taking coordinates along position first axis.

pypower.utils.exception_handler(exc_type, exc_value, exc_traceback)

Print exception with a logger.

pypower.utils.is_sequence(item)

Whether input item is a tuple or list.

pypower.utils.joint_occurences(nrealizations=128, max_occurences=None, noffset=1, default_value=0)

Return expected value of inverse counts, i.e. eq. 21 of arXiv:1912.08803.

Parameters:

  • nrealizations (int) – Number of realizations (including current realization).

  • max_occurences (int, default=None) – Maximum number of occurences (including noffset). If None, defaults to nrealizations.

  • noffset (int, default=1) – The offset added to the bitwise count, typically 0 or 1. See “zero truncated estimator” and “efficient estimator” of arXiv:1912.08803.

  • default_value (float, default=0.) – The default value of pairwise weights if the denominator is zero (defaulting to 0).

Returns:

occurences – Expected value of inverse counts.

Return type:

list

pypower.utils.mkdir(dirname)

Try to create dirname and catch OSError.

pypower.utils.pack_bitarrays(*arrays, dtype=<class 'numpy.uint64'>)

Pack bit arrays into a list of integer arrays. Inverse operation is unpack_bitarray(), i.e. unpack_bitarrays(pack_bitarrays(*arrays, dtype=dtype))``is ``arrays, whatever integer dtype is.

Parameters:

  • arrays (bool arrays) – Arrays of integers or booleans whose elements should be packed to bits.

  • dtype (string, dtype) – Type of output integer arrays.

Returns:

arrays – List of integer arrays of type dtype, representing input boolean arrays.

Return type:

list

pypower.utils.pascal_triangle(n_rows)

Compute Pascal triangle. Taken from https://stackoverflow.com/questions/24093387/pascals-triangle-for-python.

Parameters:

n_rows (int) – Number of rows in the Pascal triangle, i.e. maximum number of elements \(n\).

Returns:

triangle – List of list of binomial coefficients. The binomial coefficient \((k, n)\) is triangle[n][k].

Return type:

list

pypower.utils.popcount(*arrays)

Return number of 1 bits in each value of input array. Inspired from https://github.com/numpy/numpy/issues/16325.

pypower.utils.rebin(array, new_shape, statistic=<function sum>)

Bin an array in all axes based on the target shape, by summing or averaging. Number of output dimensions must match number of input dimensions and new axes must divide old ones.

Taken from https://stackoverflow.com/questions/8090229/resize-with-averaging-or-rebin-a-numpy-2d-array and https://nbodykit.readthedocs.io/en/latest/_modules/nbodykit/binned_statistic.html#BinnedStatistic.reindex.

Example

>>> m = np.arange(0, 100, 1).reshape((10, 10))
>>> n = rebin(m, new_shape=(5, 5), statistic=np.sum)
>>> print(n)
[[ 22 30 38 46 54]

[102 110 118 126 134] [182 190 198 206 214] [262 270 278 286 294] [342 350 358 366 374]]

pypower.utils.reformat_bitarrays(*arrays, dtype=<class 'numpy.uint64'>, copy=True)

Reformat input integer arrays into list of arrays of type dtype. If, e.g. 6 arrays of type np.uint8 are input, and dtype is np.uint32, a list of 2 arrays is returned.

Parameters:

  • arrays (integer arrays) – Arrays of integers to reformat.

  • dtype (string, dtype) – Type of output integer arrays.

  • copy (bool, default=True) – If False, avoids copy of input arrays if dtype is uint8.

Returns:

arrays – List of integer arrays of type dtype, representing input integer arrays.

Return type:

list

pypower.utils.savefig(filename, fig=None, bbox_inches='tight', pad_inches=0.1, dpi=200, **kwargs)

Save figure to filename.

Warning

Take care to close figure at the end, plt.close(fig).

Parameters:

  • filename (string) – Path where to save figure.

  • fig (matplotlib.figure.Figure, default=None) – Figure to save. Defaults to current figure.

  • kwargs (dict) – Optional arguments for matplotlib.figure.Figure.savefig().

Returns:

fig

Return type:

matplotlib.figure.Figure

pypower.utils.setup_logging(level=20, stream=<_io.TextIOWrapper name='<stdout>' mode='w' encoding='utf-8'>, filename=None, filemode='w', **kwargs)

Set up logging.

Parameters:

  • level (string, int, default=logging.INFO) – Logging level.

  • stream (_io.TextIOWrapper, default=sys.stdout) – Where to stream.

  • filename (string, default=None) – If not None stream to file name.

  • filemode (string, default='w') – Mode to open file, only used if filename is not None.

  • kwargs (dict) – Other arguments for logging.basicConfig().

pypower.utils.sky_to_cartesian(rdd, degree=True, dtype=None)

Transform distance, RA, Dec into cartesian coordinates.

Parameters:

  • rdd (array of shape (3, N), list of 3 arrays) – Right ascension, declination and distance.

  • degree (default=True) – Whether RA, Dec are in degrees (True) or radians (False).

Returns:

positions – Positions x, y, z in cartesian coordinates.

Return type:

list of 3 arrays

pypower.utils.unpack_bitarrays(*arrays)

Unpack integer arrays into a bit array. Inverse operation is pack_bitarray(), i.e. pack_bitarrays(unpack_bitarrays(*arrays), dtype=arrays.dtype)``is ``arrays.

Parameters:

arrays (integer arrays) – Arrays of integers whose elements should be unpacked to bits.

Returns:

arrays – List of boolean arrays of type np.uint8, representing input integer arrays.

Return type:

list

MPI

pypower.mpi.barrier_idle(mpicomm=mpi4py.MPI.COMM_WORLD, tag=0, sleep=1.0)

MPI barrier fonction that solves the problem that idle processes occupy 100% CPU. See: https://goo.gl/NofOO9.

pypower.mpi.domain_decompose(mpicomm, smoothing, positions1, weights1=None, positions2=None, weights2=None, boxsize=None, domain_factor=None)

Adapted from https://github.com/bccp/nbodykit/blob/master/nbodykit/algorithms/pair_counters/domain.py. Decompose positions and weights on a grid of MPI processes. Requires mpi4py and pmesh.

Parameters:

  • mpicomm (MPI communicator) – The MPI communicator.

  • smoothing (float) – The maximum Cartesian separation implied by the user’s binning.

  • positions1 (array of shape (N, 3)) – Positions in the first catalog.

  • positions2 (array of shape (N, 3), default=None) – Optionally, for cross-pair counts, positions in the second catalog. See positions1.

  • weights1 (list, array, default=None) – Optionally, weights of the first catalog.

  • weights2 (list, array, default=None) – Optionally, weights in the second catalog.

  • boxsize (array, default=None) – For periodic wrapping, the 3 side-lengths of the periodic cube.

  • domain_factor (int, default=None) – Multiply the size of the MPI mesh by this factor. If None, defaults to 2 in case boxsize is None, else (periodic wrapping) 1.

Returns:

(positions1, weights1), (positions2, weights2) – The (decomposed) set of positions and weights.

Return type:

arrays

pypower.mpi.gather(data, mpiroot=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Taken from https://github.com/bccp/nbodykit/blob/master/nbodykit/utils.py. Gather the input data array from all ranks to the specified mpiroot. This uses Gatherv, which avoids mpi4py pickling, and also avoids the 2 GB mpi4py limit for bytes using a custom datatype.

Parameters:

  • data (array_like) – The data on each rank to gather.

  • mpiroot (int, Ellipsis, default=0) – The rank number to gather the data to. If mpiroot is Ellipsis or None, broadcast the result to all ranks.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

Returns:

recvbuffer – The gathered data on mpiroot, and None otherwise.

Return type:

array_like, None

pypower.mpi.local_size(size, mpicomm=mpi4py.MPI.COMM_WORLD)

Divide global size into local (process) size.

Parameters:

  • size (int) – Global size.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

Returns:

localsize – Local size. Sum of local sizes over all processes equals global size.

Return type:

int

pypower.mpi.scatter(data, size=None, mpiroot=0, mpicomm=mpi4py.MPI.COMM_WORLD)

Taken from https://github.com/bccp/nbodykit/blob/master/nbodykit/utils.py Scatter the input data array across all ranks, assuming data is initially only on mpiroot (and None on other ranks). This uses Scatterv, which avoids mpi4py pickling, and also avoids the 2 GB mpi4py limit for bytes using a custom datatype

Parameters:

  • data (array_like or None) – On mpiroot, this gives the data to split and scatter.

  • size (int) – Length of data on current rank.

  • mpiroot (int, default=0) – The rank number that initially has the data.

  • mpicomm (MPI communicator, default=MPI.COMM_WORLD) – The MPI communicator.

Returns:

recvbuffer – The chunk of data that each rank gets.

Return type:

array_like